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ABSTRACT 

Context. The reason for the observed thinness of the solar tachocline is still not well understood. One of the explanations 
that have been proposed is that a primordial magnetic field renders the rotation uniform in the radiation zone. 
Aims. We test here the validity of this magnetic scenario through 3D numerical MHD simulations that encompass both 
the radiation zone and the convection zone. 

Methods. The numerical simulations are performed with the anelastic spherical harmonics (ASH) code. The computa- 
tional domain extends from 0.07 Rq to 0.97 Rq. 

Results. In the parameter regime we explored, a dipolar fossil field aligned with the rotation axis cannot remain confined 
in the radiation zone. When the field lines are allowed to interact with turbulent unstationary convective motions at 
the base of the convection zone, 3D effects prevent the field confinement. 

Conclusions. In agreement with previous work, we find that a dipolar fossil field, even when it is initially buried deep in- 
side the radiation zone, will spread into the convective zone. According to Ferraro's law of iso-rotation, it then imprints 
on the radiation zone the latitudinal differential rotation of the convection zone, which is not observed. 
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1. Introduction 

The discovery of the solar tachocline can be considered as 
one of the g reat achievements of helioseismology (^ee Brown 
et al.| 1989). In this layer the rotation switches from differ- 
ential (i.e., varying with latitude) in the convection zone 
to nearly uniform in the radiation zone below. Its extreme 
thinness ( less than 5 % of the solar radius) came as a sur- 
prise {see Charbonneau et allliggg ) and has still not been 
explained properly. 

'Spiegel & Zahn' (1992) made the first attempt to model 
the tachocline, and they showed that it should spread deep 
into the radiation zone owing to thermal diffusion. As a re- 
sult, the differential rotation should extend down to 0.3 Rq 
after 4.5 Gyr, contrary to what is inferred from helioseis- 
mology dSciiou et al.|1998| [Thompson et al.|2003[ ). The au- 
thors suggested one way to confine the tachocline, namely 
to counter-balance thermal diffusion by an anisotropic tur- 
bulence, which would smooth the differential rotation in 
latitude. This model was successfully simulated in 2D by 



because Miesch (2003) found in his numerical simulations 
that this turbulence is indeed anti-diffusive in the verti- 
cal direction, but the shear is reduced in latitude. Similar 
conclusions were reach ed through theoretical (analy tical) 
studies by [Kim ( 2005] ), |Leprovost fc Kim ( 2006) , and |Kim 
fc LeprovostfpQQT^ 



|Elliott| (p97D 



Gough & Mclntyre 



^_^ ( |1998D (GM98 hereafter) argued 

that such an anisotropic turbulent momentum transport is 
not able to erode a large-scale latitudinal shear, citing ex- 
amples f rom geophysical studies {see Starr_ 1968 and more 
recently [Dritschel fc Mclntyre 2008 for a comprehensive 
review). For instance, the quasi-biennal oscillation in the 
Earth's stratosphere arises from the anti-diffusive behavior 
of anisotropic turbulence. But the question is hardly settled 
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Gough fc Mclntyre proposed an alternate model, in 
which a primordial magnetic field is buried in the radia- 
tive zone and inhibits the spread of the tachocline. This is 
achieved through a balance between the confined magnetic 
field and a meridional flow pervading the base of the convec- 
tion zone. Su ch a fossil field was also invoked by Rudiger fcj 
Kitchatinov|(|1997j to explain the thinness of the tachocline 
and by Barnes et al. (1999) to interpret the lithium-7 de- 
pletion at the solar surface. Moreover, it would easily ac- 
count for the quasi uniform rotation of the radiative inte- 
rior, which the hydrodynamic model of Spiegel fc Zahn does 
not (at least in its original formulation). However, one may 
expect that such a field would spread by ohmic diffusion, 
and that it would eventually connect with the convection 
zone, thus imposing the differential rotation of that zone 
on the radiation zone below. 

In order to test the magnetic confinement scenario and 
to put constraints on the existence of such a primordial 
field, several numerical simulations of the solar radiation 
zone were perf ormed by Garaud ( 2002| ) in 2D, and by^Brunl 
fc Zahn| ( [20061 ) i^-^- paper I, hereaEer BZ06) in 3D. THese 
simulations exhibit a propagation of the differential rota- 
tion into the radiation zone (primarily in the polar region). 
The primordial magnetic field successfully connects to the 
base of the convection zone where the shear is imposed 
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and, according to Ferraro's law of iso-rotation (Ferraro 



1937), the angular velocity is transmitted along the mag- 
netic field lines from the convection zone into the radiative 
interior, much faster than it would through thermal diffu- 
sion. However, in the aforementioned calculations the base 
of the convection zone was treated as an impenetrable wall, 
i.e. radial motions could not penetrate in either direction 
(from top down or from bottom up). No provision was made 
for a meridional circulation originating in the convection 
zone that would penetrate into the radiation zone and could 
prevent the upward spread of the magnetic field. This pene- 
tration was taken i nto account in numerica l simulati ons by 
| Sule et al.l (120051), [Rudiger fc Kitchatinov] ( |2QQ7j , .Garaud 
"& Garaud (20 08|, in a model of the polar region by Wood 
\k McIntyre| ( |2QQ7| , and in simulat ions coupling the radia- 



tive zone to a convective zone by Rogers ( 2Q11[ ) to prop- 
erly represent the GM98 scenario. In most of these models, 
an axisymmetric and stationary meridional circulation was 
imposed at the top of the radiation zone, strong enough 
to bend the magnetic field lines of the primordial field. As 
could be expected, various profiles of meridional circulation 
resulted in different confinement properties, thus emphasiz- 
ing the need for a more self-consistent approach. We thus 
propose here to treat the problem by letting a genuine con- 
vective envelope dynamically generate its meridional cir- 
culation and differential rotation and to study how these 
nonlinearly generated fiows interact with a fossil field. 

Our 3D simulations treat the nonlinear coupling of the 
convection and radiation zones, and the model includes all 
physical ingredients that play a role in the tachocline con- 
finement (thermal and viscous diffusion, meridional circula- 
tion, convective penetration, solar-like stratification, mag- 
netic fields, pumping, ...). The paper is organized as follows. 
In Sect.[2]we describe the model we use to address the ques- 
tion of the magnetic confinement of the tachocline. We then 
emphasize in Sectj3]the global trend of our simulations, and 
examine in Sect. |4] the dynamics of the tachocline region. 
Discussions and conclusions are reported in Sect. [5] 



2. The model 

We use the well tested ASH code (anela stic spherica l har- 
monics, seeJClune et aT] [19991 |Miesch et al.. 2QQQ1 |Brun| 
et al. |2QQ4 ). Originally designed to model solar convec- 
tion, it has been adapted to inclu de the radiation zo ne as 

{Wn\ per- 



well (see BZ06). With this code, Brun et al. 
formed the first 3D MHD simulations of the whole sun 
from r = 0.07 i?© to r = 0.97 i?© with realistic stratifi- 
cation, thus nonlinearly coupling the radiation zone with 
the convection zone. Convective motions at the base of the 
convective envelope (around r = 0.72 i?©) penetrate into 
the radiative zone over a distance of about 0.04 i?©, ex- 
citing internal waves that propagate over the entire radia- 
tive interior. The model maintains a radiation zone in solid 
body rotation and develops a differentially rotating con- 
vect ive zone that agrees w ell with helioseismic inversions 
{see Thompson et al.|2003 ). Although approximately three 
times thicker than observed, the simulated tachocline pos- 
sesses both a latitudinal and radial shear. This tachocline 
spreads downward first because of thermal diffusion, and 
later because of viscous diffusion, as explained in SZ92. We 
do not assume here any anisotropic turbulent diffusivity to 
prevent the spread of the tachocline. 



In the present work, we start our simulation from a ma- 
ture model with a convection and a radiation zone. We then 
introduce an axisymmetric dipole-like magnetic field deep 
in the inner radiation zone. This work is meant to explore 
further the mechanisms considered in BZ06. The system is 
then evolved self-consistently throu gh the interplay of fiuid 
and magnetic field dynamics. As in Brun et aD (2011), the 



computational domain extends from VboUom = O.OTRq to 
^top = O.97i?0, and the boundary between radiation and 
convection zone is defined at r^cz ^ 0.71 5 i^p , bas ed on 
the initial entropy profile {see Figs. 1(a) and l(c)| ) that 
is assumed to setu p the background equilibrium state. We 
also define in Fig. 1(c) t he c onvective overshooting depth 
^ov ^ 0.675 i?0 {see Sect. 2.2), the penetration de pth of the 
meridional circulation rMC ^ O.GSRq {see Sect. 3.2), and 
the shear depth r shear ^ 0.58 Rq {see Sect. 2.2). 



The diffusivities (Fig. ^ were chosen to achieve a realis- 
tic rotation profile in the convective envelope (see Sect. 2.2) 
by means of an efficie nt angular momentum redistribution 



by Reynolds stresses (Brun & Toomre [2002[ Miesch et al. 
2006). In the radiation zone, they have the same values as in 
BZ06 for consistency with previous studies. We define the 
following step function to construct our diffusivity profiles: 



step{r) 



, 1 + tanh 10"^° 

2 V V ^s 

Ve^t + i^t (1 - ^e) step, 



(1) 



where Vy = 4.710^^ cm and Us = 0.1 control the radial 
localization and the thickness of the diffusivity jump. The 
parameter Vt = 8.0 10^^ cm^ s~^ controls the viscosity value 
in the convection zone and Ue = 10~^ controls the size of 
the viscosity jump. The thermal and magnetic diffusivities 
are similarly defined with the same formulation as for the 
step function, but with different hzt = 3.2 10-^^ cm^ s~^^Vt = 
1.610^^ cm^s~^ and hie = 0-25, T^g = 5.0 10~^ parameters 
to generate the diffusivity profiles plotted in Fig.|2j We used 
8i NrXNgX N^ — 1024 X 256 x 512 grid on massively parallel 
computers to compute this model. The Prandtl number v j k 
is 10~^ in the radiative zone and 0.25 in the convective 
zone; and the magnetic Prandtl number v jr] is 0.1 in the 
radiative zone and 0.5 in the convective zone. The hierarchy 
between the various diffusivities is thus respected even if 
their amplitudes are higher than in the Sun. As will be 
seen in Sect. [3J they were chosen in a way that nonlinear 
processes act efficiently in the regions of interest. 



2.1. Governing equations 

ASH solves the 3D MHD set of equations ( [2pl ). It uses the 
anelastic approximation to filter out the sound waves, and 
the LES approach (large eddy simulation) with parame- 
terization to take into account subgrid motions. The mean 
state is described by mean profiles of density p, tempera- 
ture T, pressure P, and entropy S. The reference state is 
taken from a thermally relaxed ID solar structure model 
(Brun et al. 2002) and is regularly updated. Fluctuations 
around the reference state are denoted without bars. The 
governing equations are written in electromagnetic units, in 
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Fig. 1. (a) Mean entropy gradient profile and zoom in the 
region where the entropy gradient changes sign. Entropy 
gradients are plotted in erg g~^ K~^ cm~^ . (b) Radial ve- 
locity on a spherical shell near the top of the convection 
zone (0.96 i?©). Dark colors represent the downflows, while 
bright-yellow colors denote the upflows. (c) rms velocities 
profiles at the interface between radiative and convective 
zones. The four vertical bars label the base of the convec- 
tion zone (where the entropy gradient changes sign), the 
penetration depth of meridional circulation, the overshoot- 
ing depth of convective motions, and the shear depth where 
differential rotation vanishes. 

the reference frame rotating at angular velocity flo = Qq^z- 
V • (pv) = (2) 

V-B = (3) 

p [5tv + (v ■ V) V + 2f2o X v] = - VP + pg 

+ J_(VxB) xB- V-P- [VP-pg] (4) 

pT [dtS + v-V{S + S)]=V- [KrpCpV (T + T) 



+ KopTVS + KpfVS] + ^ J2 

1 2 - 

+2piy[eijeij - - (V • v) ] + pe 

dtB = Vx(vxB)-Vx {riV x B) , 



(5) 
(6) 
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Fig. 2. Radial profile of the diffusivities used in our model 
(the values are given in cgs units) . The hz profile corresponds 
to thermal diffusivity in equation (Isl), acting on the entropy 
fluctuations. 



where v is the local velocity, B the magnetic field, k,^ v 
and T] are respectively the effective thermal diffusivity, eddy 
viscosity, and magnetic diffusivity. The radiative diffusivity 
Kj. is deduced from a ID model and adjusted in the ra- 
diative zone and in the overshooting layer to achieve an 
equilibrated radial flux balance (^ee Fig. |3|. The thermal 
diffusion coefficient k^ plays a role at the top of the convec- 
tive zone (where convective motions vanish) to ensure the 
heat transport through the surface. This term proportional 
to dS /dr is part of our subgrid scale treatment in the con- 
vection zone. The unresolved flux only acts in the upper 
part of the convection zone {see Fig. [3|, and kq is chosen 
to be small enough that this flux does not play any role 
in the radiative interior (besides eventual numerical stabi- 
lization) . pe is the rate of nuclear energy gener ation chosen 
to ha ve the correct integrated luminosity {see Brun et al. 
2011). J = (c/47r)V X B is the current density, and the 



viscous stress tensor V is defined by 



Vij = -2pv 



^^i3 ~ 3 (^ • ^) ^^i 



(7) 



The system is closed by using the linearized ideal gas law: 
p P T P S 



T -fP 



(8) 



with Cp the specific heat at constant pressure and 7 the 
adiabatic exponent. 

We chose rigid and stress-free conditions at the bound- 
ary shells for the velocity. We also imposed constant mean 
entropy gradient at the boundaries, matched the magnetic 
field to an external potential magnetic field at the top, and 
treated the bottom boundary as a perfect conductor. 



2.2. Characteristics of the background model and choice of 
the initial magnetic field 

Integrated solar models built with the ASH code that cou- 
ple a deep radiative interior t o a convective zone have been 
previously described in Brun et al. (2011); here we need 
only to stress some properties that are of particular inter- 
est for our study of the tachocline dynamics. 

The starting point of our model is the choice of thermo- 
dynamic quantity profiles that will allow us to treat the con- 
vective and radiative zones together. Using a standard solar 
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Fig. 4. (a) Azimuthally and temporally averaged rotation frequency. Yellow denotes the rotation rate of the interior, 
red is the higher rotation rate and blue the lower rotation rate, (b) Radial profile of rotation frequency at different 
latitudes, (c) Meridional circulation in the radiation zone (left) and in the convection zone (right), evaluated from the 
poloidal mass flux averaged over time and longitude. Solid contours denote counterclockwise circulation, dashed contours 
clockwise circulation. 
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Fig. 3. Radial flux balance. In the convection zone, energy 
is mainly carried by the enthalpy flux; the entropy flux 
represents the flux carried by the unresolved motions. Note 
the penetration of convective motions below the convection 
zone (represented by the dashed line). 



model calibrated to seismic data (Brun et al. 2002) com- 
puted by the CESAM code (| Morel| | 1997 p, we chose to use 



the real solar stratification in the radiation zone and take 
a constant negative initial entropy gradient in the convec- 
tion zone. These profiles are closely in agreement the guess 
profile after a Newton-Raphson solve. After setting a sta- 
ble/unstable stratification, we perturbed the background 
state assuming a Rayleigh number well above the critical 
Rayleigh number for the onset of convection. A sample of 
the convective motions realized in the model is shown in 
Fig. |l(b)] where one can recognize the characteristic banana 
pattern nea r the equator and a pat chy behavior at high lat- 
itudes (gee ||Brun fc Toomr"e||2002[ ). We then let the model 
evolve toward a mature state in two steps. First we let 



the overshooting layer and the differential rotation develop. 
Then we sped up the thermal relaxation of the system by 
adding a b ump on the radia tive flux in the overshooting 
region {see Brun et al. 2011 for more details). Thanks to 
this procedure, we quickly achieved good radial flux bal- 
ance (Fig. [3|. The energy provided by nuclear reactions is 
transported through radiation in the radiation zone, and is 
essentially transported by the enthalpy flux in the convec- 
tion zone. The overshoot region, defined here as the region 
of negative enthalpy flux, extends from rbcz = 0.715 Rq to 
rov = 0.675 i?0, leading to dov ^ 0.04 Rq ^ 0.4 Hp (where 
Hp is the local pressure scale height). In the convection 
zone, the radiative flux quickly goes to zero, while viscous 
and kinetic energy fluxes tend to oppose the outward trans- 
port of heat of the enthalpy fiux associated to the convective 
motions. Near the top of the domain, convective motions 
vanish due to our choice of impenetrable wall; thus the en- 
thalpy flux decreases, and the energy is then transported 
by the entropy flux. 

With our choice of diffusion coefficients (Fig. [2]), the 
system reached a state of fully developed convection with 
a rotation profile that agreed well with t he results of he- 
lioseismology {e.g. Thompson et al. 2003): the convection 
zone rotates faster at the equator than at the poles, a nd th e 
rotation profiles are conical at mid-latitudes (see Fig. 4(a)| ). 
We dis play t he radial profile of Q at indicated latitudes in 
Figure |4(b)| We clearly see the presence of a tachocline. 
Future work will aim at reducing the thickness of the jump 
of the diffusivities (Fig. |2]) with a more flexible radial dis- 
cretization (currently unoer development) to model a thin- 
ner hydrodynamic tachocline. Finally, the meridional cir- 
culation in the convection zone is composed by a major 
cell in each hemisphere (Fig. 4(c)[ ), and by smaller counter- 
rotating cells at the poles. 

When the hydrod ynamic model reach ed an equilibrated 
state (as reported in Brun et al. 2011), we introduced a 
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dipolar axisymmetric magnetic field buried in the radiative 
zone. It is set to vanish at the base of the tachocline, below 
the level where differential rotation starts, and we imposed 
the functional form 'B = Bq {Brer + BoGq)^ with 



Br 



t5s*, 



Be 



-A^ ■ (9) 



^(r, ^) is constant on field lines. As in BZ06, we chose 

^ = {r/Rf (r - Rb)'^ sin^ l9 for r < Rj, 
= for r>Rb 



(10) 



where R^ = 0.57 Rq is t he bounding radius of the con 



fined field. According to Gough fc Mclntyre (|1998), the 



amplitude of the magnetic field controls the scalings of the 
tachocline and the magnetopause (a thin layer of intense 
magnetic field). With our choice of parameters, these scal- 
ings would require an initial seed magnetic field of the order 
of 10^ G {i.e., S/R - A/R - IQ-^, in GM98's notations). 
Here we set 5o to 4.2 10^ G to be at equipartition of energy 
in the radiative zone between magnetic energy and total 
kinetic energy in the rotating reference frame. 

We stress here that such a mag netic field is subject to 



the high m Tayler ins tability {see jTayler 1973 Brun & 



Zahn||2QQ6[ Brun 2007) because it has no azimuthal com- 
ponent to start with. One could prefer to choose a mixed 
poloidal/toroidal magnetic field because it is the only possi- 
ble stable configuration of magnetic f ield in radiative stellar 
interiors 



Mathis 2 



(,TaylerJ973, ,Brait hwaite fc: Spruit||20Q4{ |Duez fc 
^OlOJ ). However, this has no consequences for the 



problem at hand, because a toroidal field quickly develops 
and stabilizes the magnetic configuration. 



3. General evolution 

The first evolutionary phase of our simulation is similar to 
that in BZ06, with the magnetic field diffusing through the 
tachocline. But here we follow its evolution as it enters into 
the convection zone, and we witness the back-reaction of 
convective motions on the field. The out come is the same 
in the end, as we shall see in Sect. 13.11 More details on 



tachocline dynamics are given in Sect. |3.2 



3.1. Global evolution of the magnetic field 

We started the simulation by burying the primordial field 
deep within the radiative interior and below the tachocline. 
As the simulation proceeded, the field expanded through 
the tachocline into the convection zone, where it experi- 
enced a complex evolution under the combined action of 
advection by the convective motions, shearing through the 
differential rotation and oh mic d iffusion. It finally reached 
the domain boundary (Fig. |5(b)[ ) where it connected to an 
external potential field. 

The distortion of the field lines produces magnetic 
torques that force the angular velocity ft to tend toward 
constant along the field lines of the mean poloidal field, 
thus leading to Ferraro's law of iso-rotation. This is illus- 
trated in Fig. [5J where the iso-contour Q = 414 nHz {i.e., 
the initial rotation frequency of the whole radiation zone) 
is thicker and more intense. As a result, the radiation zone 
slows down under the influence of the magnetic torques 
(Fig. 5(c)). Our choice of torque-free boundary conditions 



implies that the total angular momentum is conserved. It is 
extracted from the radiative zone and redistributed within 
the convective envelope. 

As the initial magnetic field meets the angular velocity 
shear, toroidal magnetic field is created in the tachocline re- 
gion via the l^-effect {see Moffatt (1978)), consistent with 
the magne tic la yer introduced in GM98. Near the initial 
time (Fig. |6(a)[ ) the the magnetic layer exhibits a mixed 
/ = 1,/ = 3 configuration. The axisymmetric longitudi- 
nal magnetic field then presents a / = 3 structure (Fig. 
|6(b)| ) , consistent with the action of differential rotation on 
the dipolar magnetic field through the induction equation. 
This layer of axisymmetric azimuthal magnetic field spreads 
downward, accompanying the spread of differential rotation 
(Figs. 6(c)|6(d)"] ). We also observe that the upper radial lo- 
calization of maximum azimuthal magnetic field remains 
remarkably localized at the base of our initial tachocline. 
Some azimuthal field does penetrate into the convective en- 
velope, partly because of some local l^-effect that is at work 
in this region. 

Meridional circulations exist in the radiative and con- 
vective zones and are self-consistently generated by the con- 
vective motions and local dynamics (Fig. IgI). Observe that 
the evolution of the magnetic layer greatly perturbs the 
pattern of the meridional circulations in the radiative zone. 
Comparing Fig. [6] to the hydrodynamic model (Fig. |4(c)] ), 
we observe that even if temporal averages of the meridional 
circulation produce large hemispherical cells in the convec- 
tion zone, the instantaneous pattern of this flow is multi- 
cellular both in longitude and latitude. Those 3D motions 
are unable to prevent the magnetic field from spreading into 
the convective zone, in contrast to what GM98 and Garaud] 
fc Garaud (2008) proposed in their 2D scenario. Note that 
tlogerS| ( |2011| ) confirm the absence of magnetic field confine- 
ment in recent 2D simulations, where the author coupled a 
radiative interior to a convective envelope. 

Azimuthally averaged views can mislead the reader here 
bacause from those figures one cannot deduce the real im- 
pact of the turbulent convective motions on the magnetic 
field. Figure [T] displays 3D visualizations that emphasize at 
the same time the advective and diffusive processes that act 
on the magnetic field in the convective zone. Magnetic field 
lines are twisted and sheared by convective motions, but 
still the magnetic field connects the radiative and convec- 
tive zones. One clearly observes the mostly horizontal mag- 
netic layer in the tachocline on Fig. 7(c) In that panel also 
the radiation zone starts to show the 'footprint' of the dif- 
ferential rotation of the convection zone. In order to follow 
one particular field line, we display in Fig. [8] a 3D render- 
ing time evolution using a high sampling cadence, viewed 
from the north pole. We observe that particular field lines 
(labeled A, B, and C) are strongly influenced by convective 
motions and exhibit a complex behavior, far from just sim- 
ple diffusion. The primary importance of 3D motions will 
also be discussed in Sect. 13.21 



We stress here that the magnetic field lines first pen- 
etrate into the convective zone through the equator. One 
may wonder how th is can be achieved in spite of the mag- 
netic pumping {see Tobias et al. 2001 Weiss et al. 2004) 
that should be acting in this region. Let us mention several 
reasons why the magnetic pumping may not be efficient 
enough in our simulation, and by extension perhaps in the 
Sun: 
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Fig. 5. Poloidal snapshots of azimuthal averag es of the rotation profile (colored background) and magnetic field lines 
(black). Color table is the same as in Fig. |4(a)'| 




t= 6435 d 

(d) 

Fig. 6. Poloidal snapshots of {B^^) (colored background) and instant meridional circulation (black lines). Initially, the 
magnetic field i s purely poloidal. The first figure on the left is taken just after the beginning iteration (not the same time 
as in Fig. |5(a)[ ) to demonstrate how the longitudinal magnetic field is primarily created. Red is the positive magnetic 
field, and blue the negative magnetic field. In order to display the azimuthal magnetic field in the radiation zone and the 
convection zone (where it is much weaker), a logarithmic scale was used for th e color table. The two white lines indicate 
the base of the convection zone ri)cz and the shear depth r shear {cf Fig. 1(c)). 
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Fig. 7. Three-dimensional views of the longitudinal velocity in the rotating frame. Red are the positive velocities, and 
blue the negative velocities. Picture (a) is the initial time where we added the magnetic field, and picture (b) is taken 
after two solar rotation periods. Picture (c) is taken much later in the simulation, when the magnetic layer became stable. 
Colored lines are the magnetic field lines we constructed from only a few seed points to be able to distinguish them. Dark 
colors denote lower magnetic field intensity. Note that for the third image, we only put in seed points at one longitude 
value to distinguish the field lines that explore almost the entire Lp domain in the radiation zone. 






Fig. 8. Three-dimensional views from the north pole. The left picture corresponds to Fig. |7(c) , and the two following 



ihasiz 



snapshots are taken 2 and 4 days later. We reversed the color table of the field lines to emphasize their movement. The 
green arrows indicate the movement of three field lines during this period. 



Magnetic pumping has been proven to be much less effi- 
cient when strong rotation is present and enhances hori- 
zontal mixing. This corresponds to a low Rossby number 
regime, which is the case for the Sun and our simula- 
tion. Indeed, our Rossby number Rq = uOrms/'^^ varies 
from 0.1 to 1 in the convection zone. 
It has also been proven to be much less efficient as 
the underlying region (here the tachocline) is more sta- 
bly stratified. We use here a realistic (strong) solar-like 
stratification. 

The meridional circulation tends to advect the magnetic 
field upward at the equator, opposing the effect of the 
magnetic pumping. However, our simulation is perhaps 
not turbulent enough to counterbalance this effect, even 
though Rm > 1. 

Our relatively high diffusivities may facilitate the con- 
nection between the two zones by favoring the field 
lines' expansion. Nevertheless, it is also known that dif- 



Ferreira||2009' from ideal flux frozen scenario). This im- 
plies that even though the field lines have pervaded into 
the convection zone, they may not be efficiently enough 
anchored to transport angular momentum. More details 
will be given in Sect. 13.2} 
• The nature of penetration at the base of the convection 
zone is influenced by the Peclet number. In our case, 
the Peclet number is on the order of 1 at the base of 
the convection zone, implying a regi me of oversh ooting 
rather than of penetrative convection ( ZahnllQOT ). This 
translates into a more extended region of mixing than 
what is occurring in the Sun. Hence, the motions in our 
simulation are likely to be more vigorous at rov 

For all these reasons, the magnetic pumping near the equa- 
tor in a fossil buried magnetic field scenario (Wood & 
Mclntyre 2011) is quite complex and cannot be taken for 



fusion helps the 'slipping' of the field lines (e.^., Zanni & 



granted as a mechanism to impermeably confine an inner 
magnetic field near the equat or (d etailed magnetic field evo- 
lution may be found in Sect. 3.2). 
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In order to watch how the primordial magnetic field 
manifests its presence in the convective zone, we tracked 
the evolution of Br near the surface (at r = 0.96 Rq). 
Our choice of parameters {Pm = 0.5) for this simulation 
prevents the development of a dynamo-generated magnetic 
field in the convection zone because the resulting mag- 
netic Reynolds number is less than one hundred there. The 
threshold fo r dynamo action is likely to be above that value, 
as shown by Brun et al. (2004). As a consequence, the pres- 
ence of a magnetic field at the surface of the model (as seen 
in Fig. [9]) is solely caused by the spreading and reshuffling 
of the inner magnetic field through the convective enve- 
lope. Since no dynamo action is operating in the model, 
the magnetic field decays on an ohmic time scale and will 
eventually vanish. However, we did not run the simulation 
long enough to access this decaying phase across the whole 
Sun. We also stress that the ratio ME/KE is only few 
10~^ throughout the convection zone during the late evo- 
lution of our model; it is still far too low to observe any 
infi uence of the magnetic field on the convective patterns 
{see Cattaneo et al.|20Q3 ). Note that the non-axisymmetric 
pattern of B^ on Fig. [9] is strongly correlated with Vr^ the 
vertical velocity pattern, with opposite sign of Br in the 
northern and southern hemispheres. Because the latitudi- 
nal derivative of Vr changes sign at the equator, we deduce 
that it is the shear term (B • V) v of the induction equation 
that dictates the evolution of Br in the convection zone. An 
interesting diagnostic is provided by reconstructing (with a 
potential field extrapolation) the magnetic field outside our 
simulated star, up to 2 Rq (Fig. 10). One easily recognizes 
the magnetic layer in the tachoclme and the rapidly evolv- 
ing magnetic field in the convection zone. The reconstructed 
outer magnetic field looks mainly dipolar and does not give 
any hints for its complicated inner structure. 

As mentioned above, our initial magnetic field is sub- 
ject to high m instabilities because it starts as a purely 
poloidal field. Quickly the development of a toroidal com- 
ponent stabilizes the poloidal field, resulting in a complex 
poloidal-toroidal topology. In the magnetic layer, the az- 
imuthal component of the toroidal field can locally be much 
larger than the poloidal field. S uch a toroidal field is sub- 
to low m instabilities {e.g. 
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Brun| 2007). We observe t he recu rring o: 
pattern at high latitudes on Fig. |ll(d)[ The profiles of the 
horizontal components of the magnetic field are shown in 
Figs. |1 1 (a)] and [l 1 (c) | In the two right panels of Fig. 11 , we 
subtracted the axisymmetric part of the field to render the 
instability more obvious. Note also that at the time cho- 
sen to generate the figure, the longitudinal component of 
the magnetic field is ten times stronger than its latitudinal 
component. 

3.2. Tachocline and local magnetic field evolution 

Previous studies of the GM98 scenario have shed light on 
the importance of the penetration of the meridional cir- 
culation into the tachocline (Sule et al. 2005 Garaud & 



Garaud 2008). Numerical simulations in 2D have shown 



that this large-scale fiow was able to deflect the magnetic 
fleld, and could prevent it from entering into the convec- 
tion zone. We plot in Fig. 12 the meridional circulation 
patterns realized in the bulk of the tachocline. We observe 
both the meridional circulation penetration (coming from 
above rijcz ^ 0.72 Rq), and the meridional circulation cells 




Fig. 10. Magnetic fleld lines from the center to the corona. 
For r > 0.96 Rq, a potential extrapolation is assu med. 
This snapshot is taken at the same time as Fig. 9(g) The 
two very transparent spherical shells are set at r = 0.62Rq 
and r = Rq, we display in them the contours of the non- 
axisymmetric B^ and total Br. Magenta denotes the posi- 
tive radial magnetic fleld, and yellow denotes the negative 
radial magnetic fleld. 




Fig. 11. Shell slices of the horizontal components of the 
magnetic fleld at r = 0.60 Rq in the late evolution of our 
model. Left flgures represent the full flelds, while the az- 
imuthal mean was substracted in the right flgures. Dark 
colors denote negative magnetic flelds, and bright colors 
positive magnetic flelds. 



generated in the tachocline and at the top of the radiation 
zone. In our model, the meridional circulation penetrates 
approximately 4% of the solar radius below the base of the 
convection zone. The rms strength of the meridional veloc- 
ity (shown with a logarithmic colormap) strongly decreases 
with depth in the tachocline. For instance, the meridional 
velocity loses three orders of magnitude over 0.03 solar 
radius by dropping from 18 ms~^ to 0.06 ms"^ near the 



poles at the base of the convective zone {see also Fig. 1(c) ). 
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Fig. 9. Time evolution of Vr and Br for the shell r = 0.96 Rq. The different times are exactly the same as in Fig. ^ We 
subtracted the m = component for B^. Dark colors denote negative velocities and radial magnetic field, and bright 
colors positive velocities and radial magnetic field. 



Again, because our Peclet number is on the order of 1, this 
penetration is likely overestimated. 



Meridional Circulation 



I i'l ' Jfi ' ' ^MT 'JL' 1' ' ^ 


' 11 ' EH' i 


«"!""" 1^ '^^« — "^S" -g-- — 








A J? «i W-ii k'^tvL i^¥* a 








\\rj a \\«;./J H i/''IV #I/K'^^ 


1 1 A - 'jl\^^i1 


ii"/\ *r^:' tt i// 'S^^/7:t>r:Jfi^. 


1 \wj«I^'jI^^^ 


^.s^:.'^^^ 


^^¥^_ 




psIjA^^^ci^y^-, ^ 


— V^^^^^^?SH^J=«*!^S^ir*7'-- j- \.:;^^3g|j^:^'-=1 


f;^^^^_j^ 


) \ v^ ^^^"^Mc: ^ 


i-' / J ^ 












■ "■^■^ \ 


1 [ / ^^^^^ / / / y / 


1 








# ^ 












^^ ^^N \ 


" [ i V ^k .-'■ ^■// // ^ 


r )/ - 


^V ''-''/ // ^ 


















^1^ 



latitude (deg) 

Fig. 12. Zoom of the time- averaged meridional circulation 
pattern before the introduction of the magnetic field in 
the tachocline area. The stream function is plotted in iso- 
contours, solid lines denoting clockwise circulation profiles, 
and dotted lines anti-clockwise. The background color map 
corresponds to the rms strength of the meridional circu- 
lation flow in ms~^ with a logarithmic scale. Red colors 
denote the strongest downflows. 



To evaluate the impact of the flow on the magnetic field 
evolution in the penetrative layer, we computed the mag- 
netic Reynolds number R^ = V^'^^A/rj, where V^'^^ is 
the total rms velocity, A the convective penetration depth 
under the convection zone, and r] the magne tic di ffusivity. 
We take A = ncz - Tov ^ 0.04 7?© {see Figs. 1 1(c) [ and [T2I). 
We obtained a magnetic Reynolds number on the order of 
1 to a few in this penetration zone. This means that we 
are in a regime where shear and advection are on the or- 
der of, or slightly greater than diffusion in the magnetic 
field evolution equation. Different definitions of the mag- 
netic Reynolds number may lead to different values. Indeed, 
if one retains only th e radial component of the velocity field 
(i.e., in the sense of Garaud & Garaud 2008) or only the 
meridional component of the flow, we obtain a magnetic 



Reynolds number slightly lower than one. Conversely, i f one 
uses the to tal radius of the star {i.e., in the sense of Sule 



et al.|[2QQ5 ), we obtain much higher Reynolds number val- 
ues. In the end, it is the relative amplitude and the spatial 
structure of the various terms in the induction equation 
at any given location that actually determine the evolu- 
tion of the system. The quantitative analysis of those terms 
demonstrates that we are in an advective regime, as will be 
seen in Figs. p!3p6l 

In addition, we also simulated the evolution of the same 
magnetic field in a purely diffusive case {i.e., without any 
velocity field). The evolution of the magnetic field deep in 
the radiation zone, i.e., below Tq^, is equivalent in the two 
simulations and is dominated by diffusion, its amplitude 
decreasing with the square root of time. However, the pen- 
etration of the magnetic field into the convective zone is 
faster in the full MHD case, as will be seen in Fig. [15] 

In order to demonstrate the role played by time- 
dependent 3D convective motions and nonlinearities in our 
simulation, we now turn to a detailed analysis of the evo- 
luti on of the mean poloidal and toroidal fields. Following 
Brown et al. (2010), we first evaluate the different terms 
that contribute to the creation and maintenance of the ax- 
isymmetric latitudinal magnetic field (Bq). It is essentially 
amplified by shear, transported by advection, and opposed 
by diffusion. Starting from the induction equation, one may 
write 



dtiBe) = Pfs + Pms + Pfa + Pma ^ Pc ^ P. 



MD^ 



(11) 



with Pms and Pma representing the production by the 
mean shearing and advecting flows, Pps and Pfa by fluctu- 
ating shear and advection, Pc by compression, and Pmd by 
mean ohmic diffusion. The symbol () stands for azimuthal 
average. Those six terms are 



Pps = ((B'.V)v')|, 

Pms = ((B)-V)(v)|, 

Pfa = -{(v'-V)B')|, 

Pma = -({v)-V){B)|, 



(12) 



{{vr){Bs) + {v',B',))-\np 



Pc 

Pmd = - V X (7?V X (B))| 
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These definitions are then easily transposed to the radial 
and azimuthal components of the magnetic field. 

We plot in Figs. T3p4 the mean latitudinal field ( Bp) 
and the various terms oi the right hand side of Eq. ( p!2| 
at two instants in the simulation. We first show in Fig. |13| 
the early phase of evolution of (Bq)^ and a later phase of 
evolution in Fig. 14 At radius r = 0.6 Rq (i.e., deeper 



At radius r = 0.6 Rq 
than Tov)^ we initially observe two trends. First, the mean 
advective and shear contributions effectively cancel one an- 
other, possibly because of our choice of impenetrable and 
perfect conductor lower boundary condition {{Vr) = and 
(Bq) = at r = ri)ot)- The same effect is observed for the 
(Br) component (not shown here). Second, because those 
two terms cancel each other, it is the magnetic diffusion 
that dictates the evolution of (Bq) in this deep region. This 
is not a surprise, considering our diffusivity values and the 
lack of strong motions in this deep (inner) part of the model. 




Fig. 13. Major terms of production of (Bq). The abscissa 
spans the northern hemisphere, the ordinate is a zoom in 
the tachocline. The left color bar corresponds to instanta- 
neous (Bq), while the right color bar stands for the nine 
other panels. The RHS panel is the sum of the panels e, h, 
j, and i, and represents the net evolution of the mean lati- 
tudinal magnetic field. Panels c, d, f, g, i, and j correspond 
to Pms, Pma, Pes, Pfa, Pc, and Pmd, respectively {see 
Eq. (p!2|). The panels e and h correspond to the sum of 
panelsc and d, and f and g, respectively. Color levels are 
not linear in order to see at the same time the contributions 
of the terms in the different regions. Red colors denote the 
positive contribution to the magnetic field and blue colors 
the negative contributions. The two horizontal lines in each 
panel represent rijcz (dash) and vmc (dot-dash), as in Fig. 



On the other hand, the patterns are much more com- 
plicated in the penetration layer (between Vov = 0.675 Rq 
and ri)cz ^ 0.72 Rq). Mean advection and shear are impor- 
tant, their sum (panel e) mainly acts at latitudes higher 
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Fig. 14. Major terms of production of (Bq). The panels 
represent the same as in Fig. [13) but at a later time. 



than 40°. Diffusion (panel j) acts everywhere in the pen- 
etration layer, but its action does not dominate the other 
contributions (panel b). The non-axisymmetric contribu- 
tions (panels f, g, and h) are very important at latitudes 
lower than 60 in the penetration layer. Therefore, 3D mo- 
tions are obviously responsible in our model for the fail- 
ure of the magnetic field confinement below the convection 
zone at low latitude. This result contrasts with previous 2D 
studies, where these motions were not taken into account. 



During the late evolution of our model (Fig. 14), 
the complex balance between axisymmetric motions, non- 
axisymmetric motions and diffusion is still operating in the 
penetration layer. We observe (panel a) that the magnetic 
field has completely pervaded the tachocline and extends 
into the convection zone. Again, the relative importance of 
the three panels e, h, and j indicates that diffusion is not 
the dominant process transporting the field. Indeed, we find 
again that the non-axisymmetric components significantly 
act at low latitude in this layer. We also note here that the 
levels of the color table were modified between figures [13] 
making the diffusive term (which is still operating 
. r = 0.6 Rq) harder to be seen in panel j. Finally, we 



14 



and 

arouna 

note that magnetic pumping is proven inefficient to confine 

the magnetic field at both times. 

In order to estimate the pumping, we plot in Fig. [l5 
the evolution of the amplitude of the latitudinal magnetic 
field at r = 0.86 Rq (in the convection zone) after its intro- 
duction. As previously seen, the amplitude rises. We plot 
on the same figure the evolution of the absolute value of 
the latitudinal component of the magnetic field both in the 
downflows (dot-dash) and in the upflows (triple dot-dash). 
We observe that the amplitude rises slower in the down- 
flows than in the upflows. In addition, the temporal lag 
between the up- and downflow regions increases with time. 
This agrees with the action of magnetic pumping by down- 
flows. However, it is a minor effect in the global evolution 
of the magnetic field. 
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Fig. 15. Rms and absolute value of the latitudinal mag- 
netic field amplitude evolution at r = 0.86 Rq. The dotted 
and dashed lines correspond to the absolute value of the 
latitudinal component of the magnetic field in the down- 
fiows and in the upflows. We selected the downflows part 
by setting the magnetic field to zero at the points corre- 
sponding to a positive Vr on the r = 0.86 Rq spherical 
shell. We averaged the remaining values over the shell to 
obtain Bg in downflows, and did the opposite for the up- 
flows selection. The diamond line represents the evolution 
of the rms latitudinal magnetic field in the purely diffusive 
(test) case. 



We also plot in Fig. 15 the evolution of the rms lati- 
tudinal component of the magnetic field in the convection 
zone in both the nonlinear simulation discussed in the paper 
(solid line) and in another purely diffusive test case (dia- 
monds). Once more, this stresses that the field evolution is 
different in our simulation than in the purely diffusive test 
case. The larger amplitude of Bq in our convection simu- 
lation clearly demonstrates the role of advection in pulling 
the field inside the convection zone. The trends in Fig. 15 
are similar between r = r^^ and rtop^ i-e.^ in the whole 
convection zone and even in the penetration layer. 

We may now also examine two different processes in 
the evolution of the azimuthal magnetic field: the creation 
and sustainment of the magnetic layer, and the expansion 
of the magnetic field into the convection zone. {B^^p) is 
firstly created underneath the convection zone through 
the interaction between the dipolar magnetic field and the 
differential rotation around r = 0.62 Rq. Diffusion slightly 
counterbalances the Q effect, and advection does not act 
much in this layer. In order to illustrate these processes, 
we plot in Fig. 16 (B^) and the right hand side terms 



from Eq. ( |12| (as m Fig. 13, but taken in the cp direction) 
during the early evolution of our model. Diffusion (panel 
j) opposes the creation of magnetic field by mean shear 
(panel c) in the magnetic layer (around 0.6 Rq). Although 
the magnetic layer seems to remain radially localized, 
the magnetic field also exists in the convective zone. 
It has a much smaller amplitude than in the magnetic 
layer, but it is not negligible. We observe that (B^p) is 
preferentially created through the two terms of shear at 
latitudes below 40°. Notice that the non-axisymmetric 
velocity and magnetic field (panels f, g, and h) contribute 
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Fig. 16. Major terms of production of {B^p) in the 
tachocline. The various contributions are taken at the same 
time as Fig. [ISJ but using expressions corresponding to 



equally with the axisymmetric terms to produce some 
magnetic field in the penetration layer (between the two 
horizontal lines). Diffusion generally tends to oppose the 
mean shear contribution Pms at all latitudes at the base of 
the convection zone, but does not act much in comparison 
with the other contributions in panels e and h. The mean 
advection opposes the creation of negative {B^p) above 
latitude 40°, which explains why the significant creation 
of coherent mean azimuthal magnetic field is localized at 
lower latitudes. 

To summarize, we showed that the magnetic field evo- 
lution through the tachocline is caused by a complex inter- 
play between shear and advection mechanisms. The diffu- 
sive process does not dictate the evolution of the magnetic 
field at the base of the convection zone, neither at the equa- 
tor nor in the polar regions. Magnetic pumping is present, 
but inefficient. As a result, the magnetic field penetrates 
into the convection zone and does not remain confined. We 
stress here that the nonlinear term involving Maxwell and 
Reynolds stresses is a key and fundamental feature of our 
3D simulation. 



4. Balances in the tachocline 

In order to understand the global rotation profile realized in 
our simulation {see Sect.|3|, we study the MHD meridional 
force balance and the angular momentum transport equa- 
tion in the following subsections. Even if the magne tic fi eld 
does not modify the meridional force balance (Sect. [XT] ), it 
is proven to be responsible for the outward transport of an- 
gula r momentum in the mid- to high-latitude regions (Sect. 
4.2). 
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4.1. MHD meridional force balance 



Starting from the vorticity equation (Eq. (13)), one can 
derive the fuh meridional force balance in the general mag- 
netic case (e.^., |Fearn||1998[ |Brun]|2QQ5D . 



dfU? 



-V X {iVa X V) 



1 



-Vp X VP- V X 



P9 



-5-Vx ( -(VxB) xB 
47r \p 



Vx ( iv-P) ,(13) 



with cJa = V X V + 2fto the absolute vorticity and u; = 
V X V the vorticity in the rotating frame. In a statistically 
stationary state, equation (13) can be rewritten by taking 
the temporal and azimuthal averages of the longitudinal 
component: 
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In Eq. (14) we underlined different terms: 



stretching represents the stretching of the vorticity by 
velocity gradients; 

advection represents the advection of the vorticity by 
the flow; 

compressibility represents the flow compressibility reac- 
tion on the vorticity; 



^-^ is the baroclinic term, characteristic of non- 
aligned density and pressure gradients; 
Vdc~ ^ ae ^^ P^^^ ^^ ^^^ baroclinic term but arises from 
departures from adiabatic stratification; 
viscous stresses represent the diffusion of vorticity 
caused by viscous effects; 

magnetic torques 1 and 2 represent the 'shear' and 
'transport' of the magnetic field by the current; 
compressibility 2 represents the component of the curl 
of the magnetic torque owing to compressibility of the 
flow. 



The meridional force balance equation (14) helps under- 
standing which physical effects are responsible for breaking 
the Taylor-Proudman constraint of cylindrical differential 
rotation in the conv ection zone (a chieved for example with 
barochn ic flows, see [Miesch et al.||2QQ6| [Balbus et ar|[2QQ9] 
[Brun et al.|2QlQ ). Indeed, considering a small Rossby num- 
ber and neglecting stratiflcation, Reynolds, Maxwell, and 
viscous stresses, Eq. (14) reduces to 



d{v^) ^ g d{S) 
dz 2rtorCp dO 



(15) 



which is the origina l thermal wind equation {e.g. Pedlosky 
1987; Durney 1999). The baroclinic term induces motions 
that break the cylindrica l rotation proflle. In the model of 



Wood & Mclntyre (2011), this balance is considered in the 
polar regions. Compositional gradients and magnetic effects 



are added to equation (15). This balance is of primary im- 



portance in their magnetic conflnement layer scaling. We 
therefore present below the full meridional force balance of 
our model to evaluate the role played by all mechanisms in 
the tachocline, except compositional effects. 

Fig. 17 displays the main contributing terms from Eq. 



( 14 ) in our simulation. The magnetic terms were increased 
(iBy a factor of 500) to make them more visible. Surprisingly, 
they are too low to contribute to the meridio nal force bal - 



ance, as such there is n o 'magnetic wind' ( gee [Fearn||1998 ). 
As already noticed in Brun et al. (2011), the baroclinic 



term dominates the balance, especially under 0.65 i?©, ie., 
below the penetration. At the base of the convection zone 
(0.72 i?0), all terms play a signiflcant role. Stretching, ad- 
vection, and viscous stresses compete against each other at 
low latitude, letting the baroclinic term dominate again. 
Near the pole, the balance is more complex, and viscous 
stresses (helped a little by stretching and advection) coun- 
teract the positive/negative pattern of the baroclinic term. 
The strict thermal wind balance expressed in Eq. (15) is 



thus not observed in the polar region. The magnetic terms 
are located in the magnetic layer where B^^ is generated 
{see Sect.pJ, and are negligible. Indeed, our magnetic fleld 
is far from being force-free, and the Lorentz force is non- 
zero in the three components of the momentum equation. 
The poloidal velocity is mainly driven by advection of con- 
vective motions and by baroclinic terms in the radiation 
zone. As a result, the poloidal magnetic fleld does not act 
much on the poloidal velocity. Note that this does not imply 
the opposite. We indeed observed {see Sect. 3.2) that the 
poloidal magnetic fleld evolution is dominated by advection 
and shear fro m convective motion s . We d id not recover the 
balance from Wood fc Mclntyr'e] (2011), although we al- 



ready have a sufficiently strong magnetic fleld to transport 
angular momentum, as seen in Sect. |4.2[ Gyroscopic pump- 
ing is not operating efflciently near the polar region. There 
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Fig. 17. Zoom around the tachocline of the different terms 
discussed in Eq. (14), averaged over longitude and two solar 
rotation periods during the late evolution of the model. 
Magnetic terms were increased by a factor 500 to make 
their structure more evident. 



is a small downwelling, b ut it is not s trong enough to modify 
the local dynamics (gee|Brun et al.||2011_for more details). 



4.2. Angular momentum balance 

Identifying processes that redistribute angular momentum 
in the bulk of the tachocline gives us clues to under- 
stand how this layer evolves. Because we choose stress- free 
and torque free magnetic boundary conditions, no exter- 
nal torque is applied to the system. As a consequence, 
the angular momentum is conserved. Following previous 
studies (e.^., Brun et al.| |2004), we examine the contri- 
bution of different terms in the balance of angular mo- 
mentum. Averaging the (f component of the momentum 
equation Q over (f and multiplying it by rsin^, we obtain 
the following equation for the specific angular momentum 
£ = rsin6> {ftorsinO + {vcp)): 



dtipC) = -V • (F^^ + F^^ + F^^ + F^^ + F^^), (16) 



hydro 



MHD 



where the different terms correspond to contributions from 
meridional circulation, reynolds stress, viscous diffusion. 



;^MC _ 
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sm 
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r sin^ 
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r sin^ 

47r 



{B^)(Bm) 



(17) 
(18) 



(19) 

(20) 
(21) 



where the subscript .m designates the meridional compo- 
nent of V and B. In the previous equations, we decomposed 
the velocity and the magnetic field into an azimuthally av- 
eraged part (.) and (/^-dependent part (with a prime). The 
different contributions can be separated between radial (J> 
along e^) and latitudinal {J'e along e^i) contributions. We 
then compute a radial flux of angular momentum defined 

by 



Mr) 






Fr{r,e)r'^sme dO, 



(22) 



where (6>i, O2) maybe chosen to study a particular region of 
our simulation. 

As emphasized in Sect. |3.1[ magnetic field lines seem 
to be advected and diffused first outward and then pole- 
ward. In Fig. [18] we plot the different components of the 
radial angular momentum flux in the radiative zone and 
the tachocline at different times, summed over the polar 
region and at mid latitudes. Initially, the magnetic field 
is axisymmetric so Maxwell stresses do not contribute to 
the angular momentum balance. Because we chose a purely 
meridional magnetic field, there is no large-scale Maxwell 
torque. The radial angular momentum flux is then essen- 
tially carried by the meridional circulation in the radiative 
zone and is very weak (few percents of its value in the 
convective zone). Note, however, that internal waves are 
present in our simulation and m ay contribute to t he trans- 
port of angular momentum (see Brun et al.||2011 for more 
details). In the tachocline, meridional circulation, Reynolds 
stress and viscous diffusion transport the angular momen- 
tum. 

While the magnetic field more and more pervades the 
domain, a magnetic torque develops in the magnetic layer 
(along with B^ generation, see Figs. [5][6| at the top of 
the radiative zone both at high and miolatitudes. Angular 
momentum is then extracted by the Maxwell torque from 
the radiative zone (that locally rotates faster) into the 
tachocline (that locally rotates slower), slowing down the 
radiative interior. This torque tends to align the Vt con- 
tours with the poloidal magnetic field lines, thus bringing 
the system closer and closer to a state known as Ferraro's 
law of isorotation. We remark in addition that the Maxwell 
torque operates initially in a wider part of the radiative 
zone in the polar region than at mid latitudes. Because the 
magnetic field primarily connects at the poles, it is natural 
that the transport of angular momentum by the magnetic 
field starts earlier in this region. One may also observe that 
this transport is more intense at mid-latitude, acknowledg- 
ing that the magnetic field lines involve a magnetic field 
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Fig. 18. Angular momentum redistribution in the radiative zone and the tachocline region. Panels (a) are summed over 
Oi = 60 and O2 = 90 (north pole), and panels (b) are summed over Oi = 30 and O2 = 60 (mid latitudes). The different 
fluxes are averaged approxima tely over six rotation periods. From left to right, the three plots correspond to the instants 
used to make Figs. 6(a), 6(c) and|6(d)] Our magnetic field does not initially change the angular momentum balance. 



with much higher strength there. We stress again that the 
angular momentum is conserved in our simulation, mean- 
ing that the angular momentum lost in the radiative zone 
is redistributed in the convective zone. 

We also observe that the magnetic torque applies essen- 
tially in the radiation zone and becomes negligible at the 
base of the convection zone. Because we are not yet in iso- 



rotation {see Fig. 5(d)| , this may only be interpreted by 
ak 



the fact that B^ is weak in this region. 

Because the first magnetic field lines that transport an- 
gular momentum are localized at mid- and high-latitudes, 
we observe here only a slow-down of the radiative interior. 
If we had continued this simulation, magnetic field lines 
would have connected the radiative zone to the convection 
zone near the equator. This would imply a speed-up of the 
radiative interior, consistent with the results from [Brun fc 



|Zahn| ([2006). However, we do not show this later evolution 
because it does not add any useful information. 

Viscosity seems to act substantially in the bulk of the 
tachocline, augmenting the angular momentum transport 
locally. This viscous transport is caused by the radial shear 
and the relatively high viscosity value used to meet numer- 
ical requirements. As a result, diffusive processes are over- 
evaluated in our simulation. Angular momentum transport 
like this would not contribute in the Sun, given the low 
viscosity of the solar plasma. In our case, both thermal 
and viscous diffusion processes are acting in our simulation 
thanks to our choice of Prandtl number (10~^) in the radia- 
tive zone. However, BZ06 showed that the simulated viscous 



spreading of the tachocline does not make much qu alitative 
difference with the thermal spreading introduced by [Spiegel 
[fc Zahn| ( 119921 ). 

The typical viscous diffusive time scale in our simula- 
tion is given by r^ = R^ /v ~ 2.2 10^^ s at the top of the 
radiative zone. As the magnetic field grows, the Alfven time 
ta = R^Aiip/ B diminishes alike and magnetic effects be- 
come more important than diffusion. The angular momen- 
tum balance analysis emphasizes two phases in our simu- 
lation: a diffusive transport of angular momentum and the 
action of large-scale Maxwell torque. This confirms that our 
magnetic field does not prevent the spread of the tachocline, 
but on the contrary transports angular momentum out- 
wards and speeds the process up. 

5. Discussion and conclusion 

The existence of a magnetic field in the solar radiation 
zone has been the subject of a longstanding debate. Such 
a field, for example of fossil origin, could easily account 
for the uniform rotation of the radiative interior. It could 
also prevent the inward spread of the tachocline, through 
thermal diffusion, as was proposed by |Gough fc Mclntyre| 
(1998). However, it has been argued that such a field would 
diffuse into the convection zone and imprint the latitude- 
dependent rotation of that region on the radiation zone be- 
low (which is not observed by helioseismology inversions). 
This obje ction ha s indeed bee n confirmed by t he cal - 
culations of |Garaud| ( [2002, ) and [Brun fc Zahn] ( |2006 ). 
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However, their simulations treated the top of the radiation 
zone as a non-penetrative boundary, which does not allow 
for a meridional flow that would originate in the convection 
zone and could oppose the diffusion of the magnetic field. 
This r estriction was lif t ed in the 2D simulations per - 
formed bylSule et al.|(|2QQ5D ;|Rudige r fc KitchatiiiQvl (|2QQ7D; 
Garaud & Garaud ( 2008 ) , who showed that this merid- 
ional circulation could indeed prevent the magnetic field 
from diffusing into the convection zone. Very recent work 
in 2D (Rogersl2011) has studied the role of a fossil magnetic 



field in a simulation coupling a radiative and a convective 
zone. It confirms our result that the field is not confined 
at the base of the convection zone. However, the author 
does not find any magnetic transport of angular momen- 
tum, possibly because of a force-free state for the magnetic 
field. These 2D simulations assumed large-scale axisymmet- 
ric flows, unlike the three-dimensional motions of the real 
convection zone. For this reason we decided to treat the 
problem using a 3D code, with the best possible represen- 
tation of the turbulent convective motions. 

We showed by considering time-dependent 3D motions 
that the interior magnetic field does not stay confined in 
the radiative interior, but expands into the convection zone. 
Non-axisymmetric convective motions interact much differ- 
ently with a dipolar-type magnetic field, as would a mono- 
lithic, single cell large-scale meridional circulation. 

That our magnetic field is able to expand from the ra- 
diative interior into the convection zone results in an ef- 
ficient outward transport of angular momentum through 
the tachocline, imposing Ferraro's law of iso-rotation. It is 
likely a direct consequence of the non force-free configura- 
tion of our magnetic field that yields this efficient transport 
through the large-scale magnetic torque. 

We are aware that our simulations are far from repre- 
senting the real Sun, because of the enhanced diffusivities 
we had to use to meet numerical requirements. The over- 
shooting convective motions are overestimated owing to our 
enhanced thermal diffusivity Hi. The high solar Peclet num- 
ber (lower hi) regime would make the flow amplitude drop 
even faster at the base of the convection zone, which would 
result in an even weaker meridional circulation velocity. 
This would not help the confinement of the magnetic field. 
Although our enhanced magnetic diffusivity is quite large, 
we demonstrated {see Sect. |3| that our choice of param- 
eters ensured that we are in the correct regime for mag- 
netic field confinement, with advection dominating diffu- 
sion. Nevertheless, more efforts will be made in future 
work to increase the spatial resolution and to improve the 
subgrid-scales representation. 

The initial magnetic topology may also be questioned. The 
gradient of our primordial dipolar field possesses a maxi- 
mum at the equator. Because the field is initialized in the 
radiation zone where its evolution is diffusive, the magnetic 
field initially evolves faster at its maximum gradient loca- 
tion, i.e., at the equator. As a result, it first encounters 
the complex motions coming from the convection zone at 
the equator. Other simulations conducted with a modified 
topology to change the location of the maximum gradient 
(and thus the latitude where the magnetic field primar- 
ily enters the tachocline) led to similar results as those re- 
ported here. 

Given the extension of our tachocline, the differential 
rotation is already initially present below the base of the 
convection zone before we introduce the magnetic field. A 



confined magnetic field would naturally erode this angular 
velocity gradient, though we do not know much about the 
real initial conditions. In our simulation the magnetic field 
opens into the convection zone, consequently there is no 
chance to obtain a solid rotator in the radiative zone be- 
cause there is no confinement. In fact, our magnetic field 
lines primarily open at the equator. We point out that this 
is a major difficulty with the GM98 scenario: even if mag- 
netic pumping can be a way to slow down the penetration 
of the magnetic field into the convective zone, it will not 
prevent it from connecting the two zones, and/or prevent 
the magnetic field lines from opening there. Furthermore, 
the meridional circulation is not a laminar uni-cellular flow 
in our 3D simulation. It has a complex instantaneous multi- 
cellular pattern, which affects the efficiency of a polar con- 
finement. 

We therefore conclude that a dipolar axisymmetric fos- 
sil magnetic field cannot prevent the spread of the solar 
tachocline. We intend to explore other mechanisms and dif- 
ferent magnetic field topologies. Our 3D model will allow 
us to test non-axisymmetric configurations, e.g., an oblique 
and confined dipole. 
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